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ABSTRACT 



A recurring problem in solid state physics is that of 
determining the electrostatic field or potential at a point 
within the crystal lattice arising from the molecules of the 
lattice themselves. A typical problem is that of finding 
the electrostatic binding potential of an ionic lattice, 
assuming point changes. 

Nijboer, de Wette and others have worked toward 
perfecting a method for calculating the electric field 
arising from a lattice of dipoles. This work culminated 
in a paper by de Wette and Schacher detailing a procedure 
which ensures a rapidly converging summation of the three 
dimensional lattice sums involved. 

The present paper is an extension of the latter, 
presenting a computer program for evaluating dipole fields 
in a crystal lattice using the equations obtained by de Wette 



and Schacher. 
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V 



BACKGROUND 



A recurring problem in solid state physics is that of 

determining the electrostatic field or potential at a point 

within the crystal lattice arising from the molecules of the 

lattice themselves. A typical problem is that of finding 

the electrostatic binding potential of an ionic lattice, 

assuming point changes. 

1 2 

Nijboer, de Wette * and others have worked toward 

perfecting a method for calculating the electric field 

arising from a lattice of dipoles. This work culminated 

-a 

in a paper by de Wette and Schacher^ detailing a procedure 
which ensures a rapidly converging summation of the three 
dimensional lattice sums involved. 

The present paper is an extension of the latter, 
presenting a computer program for evaluating dipole fields 
in a crystal lattice using the equations obtained by de Wette 
and Schacher. 
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B.R.A. Nijboer and F.W. de Wette, Physica 14, 422 
(1958) 

^F.W. de Wette, Phys. Rev. 123 » 103 (1961) 
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F.W. de Wette and G.E. Schacher, Phys. Rev. 137, 
A?8 (1965) 



THEORY* 4 - 



Virtually all problems concerning crystal self fields 
reduce ultimately to evaluating a lattice sum of the 
general form: 



OO 



S (X,J^JL) =H fiXrOQj , \L-ba x ^-COg) 

c»,b,c=-oo 



where x, y, and z are coordinates of the point at which 
the sum is to be evaluated; f ( ) represents the contri- 

bution to the sum from a particular lattice point; and 
aa-j_, bag, ca^ are the coordinates of the lattice points 
whose contribution is being considered (a-^, a-,, and a^ 
being lattice basis vectors). 

Finding the electric field due to a lattice of dipoles 
reduces to evaluating such a sum, where the summation 
is, mathematically, conditionally convergent. That is, the 
value of the sum depends on the order in which the sum- 
mation is carried out ( analagous to the fact that the 
internal field of a crystal depends on the shape of the 
crystal). Thus, any procedure for evaluating dipole sums 
must, at some point, make an assumption of a particular 
crystal shape. 

^For a complete description see (2) 
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In the procedure used here these conditionally 
convergent triple sums over the real lattice are converted 
into absolutely (rapidly) convergent double sums over the 
recipricol lattice. Converting from a triple to a double 
sum forces an assumption of a slab-shaped crystal with 
the short dimension in the z-direction (perpendicular to 
the plane described' by the lattice vectors a_^ and a^ • 

See figure 1 for relations between x, y, z coordinates 
and the lattice basis vectors) . 

To achieve the conversion to an absolutely convergent 
double sum, the original sum is first Fourier transformed, 
making use of the property of lattice sums that the sum- 
mation of a function over a lattice in real space is equal 
to the sum of its Fourier transform over the recipricol 
lattice (in Fourier space). In Fourier space it is then 
possible to carry out the z-summation analytically, leaving 
an absolutely convergent double sum in recipricol space 
(with the assumption of a slab-shaped form). 

A. distinction is made between the cases where the 
evaluation point does or does not lie in an x-y plane of 
dipoles (.z = ca^ or z ca^ , c an integer). In the first 
case the function to be summed has a simple pole (with 
respect to the z-summation) at z = 0. In this event an 
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auxilliary convergence function is used and the sum is 
broken up into two sums, one in real space over the x-y 
plane containing the evaluation point and one in recipricol 
space over the rest of the lattice. In the second case 
the function is everywhere regular and is Fourier trans- 
formed as it stands. 

In either event the method is directly applicable only 
to lattices in which all dipoles are of the same magnitude 
and orientation. To evaluate the field of a lattice contain- 
ing more than one orientation and/or magnitude of dipoles 
it is necessary to divide the lattice into sublattices, each 
containing only one kind of dipole, evaluate the field for 
each sublattice, and combine the fields by superposition. 

The summation procedure, applied to a sublattice, 
produces nine numbers S-. s i,j = x,y,z ; where Su is the 
i-component of the electric field (in Coulombs) due to a 
lattice of unit dipoles (one electronic charge and length 
1A) oriented in the j-direction. 

While program inputs are in terms of the normal 
crystallographic quantities (a^, a-,, a^ , c*, y, S ) the 
calculations and outputs are in terms of a cartesian 
reference frame related to the crystal basis vectors as 
detailed in figures 1-3* 



4 




Figure 1: Unit Cell: 1 dipole per unit cell volume, 

Basis vectors a^, ag, a^. Crystal angles s <X 9 ^.£U,a t ; 
/3-Z^^, ; %*X-£L,Q+ ; 

x axis parallel to a^, y axis in a^ — ag plane perpendicular 
to x axis, z axis perpendicular to x-y plane. 
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Figure 2: Point at which field is to be evaluated: 

located with respect to lattice point by fractions of 
basis vectors Q +jxSk J 0 ^ 
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Figure 3s Decomposition of arbitrary dipole: 
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Then given a lattice of arbitrary dipoles £, as in figure 
3, we have: 

p x = p sin 9 cos 0 
Py = p sin O sin (J> 
p z = p cos e 

so that the field at some point j ( j , jg, j^) is given by: 

E x (j) — Px®xx^'^ "*"^y^xy^*^ "*"Pz®xz^ 

E y (j) = P X S yx (j)+PyS yy (j)+P 2 Sy Z (j) 

E z0) = Px S 2 x(j) + Py S zy(i> + pz S zz(j) 

The nine S-- are not all independent of course. They 

satisfy 

S ij = S ji 

s xx + s yy + s zz = 0 

for all lattices. Other special relations (such as S.j = 0, 
i # j, for cubic lattices) are also obtained for lattices 
with a high order of symmetry. 

In a lattice containing more than one kind of dipole 
two special considerations arise in decomposing the lattice 
into sublattices. First and most obviously, one must 
choose the j's in each sublattice, determining the point at 
which the fields are to be evaluated, such that they refer 
to a common point in the composite lattice. The second 
consideration refers back to the shape dependence of the 
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sum. In computing the sums over each sublattice, the 
shape dependence has introduced a depolarization factor 

O Q 

L = ^/(vol of unit cell in A- 3 ) into the S sum. Then 

z z 

to superpose the fields from several sublattices with 
various z-directions (referred to the composite lattice) 
it is necessary to modify all the Sjj to bring all sublattices 
to a common x-y-z orientation by manipulating the 
depolarization factors for each sublattice. A similar 
consideration arises in comparing fields obtained by this 
method and those obtained by some other method which 
assumes a different crystal shape or with experimental 
results on a specific crystal of known shape. Kittel 
considers the general problem of depolarization factors 
in some detail and tabulates the depolarization factors 
for most commonly used shapes (spherical, cubic, slab, 
etc.) An example of such considerations for a benzene 
model is contained in Appendix D. 
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C. Kittel, Introduction to Solid State Physics, 

2nd edition, 158, (dJohn Wiley and Sons, Inc., New York, 
1956) 
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PROGRAM COMPOSITION 

The actual translation of the method of Schacher 
and de W«tte is relatively straightforward although quite 
lengthy. 

The program is written in the FORTRAN 60 language 
for the CDC 1604 computer and is compatable with most 
FORTRAN monitors. The program is divided into a main 
program, containing all input-output routines and a segment 
to choose a subroutine, and eight subroutines to do the 
actual summation. 

Simplified flow charts of the entire program and 
various segments are shown in figures 4-7 • 



Appendix A contains the program listing and details 
of timing, storage, etc. 
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Figure h'* Flow Main Program - LATSUM 



if _error. 
i in input 



e 



utput error 
code 







SI 
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Figure 5: Common Flow-Subroutine SUM1, SUM2, SUM3 



Main program uses 



SUM1 if XC3^0, 
SUM2 XJ3/0, 
SUM3 XvJ3^0 , 



xdel/o 

XDEL=0 , XJ1+XJ2/0 
XDEL=0,XJ1+XJ2=0 
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Figure 6: Common Flow - Subroutine SUM2+, SUM5, GAM1 

Main program uses SUM4 &■ SUM6 if XJ3=0, XDEL^O 

SUM5 & SUM6 if XJ3=0, XDEL=0 



SUM4 or SUM5 
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Figure 7: Flow Subroutine SUM6, GAM2 

Used whenever X«J3=0 



SUM6 
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As can be seen from the flow charts, a deck can be 
made up for a particular lattice sum by appending to the 
main program either subroutines SUM1, SUM2, and SUM3 
or subroutines SUM1+, SUM5, SUM6, GAM1 and GAM2 
depending on whether is not zero or zero respectively. 
Further, if sums are to be computed for several sets of 
data the main program deck can be modified by looping 
from the input segment to the END card. 

Inputs to the program are given in Table 1. 

Output is scattered throughout the main program. 

A check output of the input parameters always appears. 

If there is an error in the input parameters the following 
is printed and control returned to the monitor: 

PROGRAM STOPPED BY INPUT ERROR 
DESIGNATED BY K= 

The individual error codes are: 



K ERROR 



1 


A1 — 0 . 




2 


A2 — 0 . 




4 


A3^ 0. 




10 


XvJ 1 < 0 . or 


XJ1 >1. 


20 


X J2 < 0 . or 


XJ2 > 1. 


40 


X«J3<0. or 


XJ3> 1. 


100 


XALPH 0 , 


. or XALPH > ' 


200 


XBET^ 0 . 


or XBET > 90. 


400 


XGAM ^ 0 . 


or XGAM ^180 


1000 


XDEL < 0. 


or XDEL > 90. 


2000 


N — 0 
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Table 1: Program Inputs 



LABEL 


UNITS 


DESCRIPTION 


ALLOWED 

LIMITS 


A1 

A2 

A3 


A 


Length of basis vector aj 


aj> 0 


XJ1 

X<J2 

X«J3 


Dimen 

sion 

less 


jj locating pt at which field 
is to be evaluated in terms of 
ai (see figure 2) 


3 ~ Ji ~ 1 


XALPH 

XBET 

XGAM 


Degrees 

and 

decimals 


Crystal angles (see fig 1) 


3 < XALPH^ 90 
KXBET^90 
D < XGAM< 180 


XDEL 


degrees 

and 

decimals 


Inclination of A3 axis from 
perpendicular .(see figure 1) 
Enter value if known else 
enter XDEL=90 and program 
will compute 


) ^ XDEL ^90 


AEX 

REX 

AEY 

REY 

AES 

RES 


A=>/?"' 

R=£ 

dimen 

sion 

less 


Error limits for convergence 
tests. Summation will contin- 
ue until the next term to be 
summed is less than 
(x) (x) 

AE( y) +RE( y) times absolute 
( z ) ( z ) 

value of sum at that time 

AEX&REX control Sj x &-S x j 

AEY&REY control Sj &-S*. 

AE2&RE2 control S,,&S , 

iz zl 

\ fuller discussion is given in 
Appendix B 


NONE 


N 


dimen 

sion 

less 


Summation cutoff. Stops sum- 
nation after covering 2N by 2N 
;erms in x-y plane if converg- 
ence t.g.st.?. not rpep fref.g.Eg , then 


N > 0 
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Multiple errors result in a K value which is the sum of the 
individual K- values involved. 

A second error output (and transfer back of control) 
occurs if N is so small that it is impossible to check the 
sums for convergence. This subject is covered more 
completely in appendix B. The printout is self-explanatory: 

PROGRAM STOPPED BECAUSE N TOO SMALL. 

N MUST BE AT LEAST + 1 TO ENSURE A 

GOOD TEST FOR CONVERGENCE. 

These are the only error outputs contained in the 
program other than those that might arise from library 
routines. There are no error stops, pauses, dumps, 
sense switches or sense lights used. 

Normal output (presuming no input errors) includes 
the nine sums Sjj, displayed as a labeled array. Following 
this is one of the statements: 

(SUM1) 

SUMMING WAS DONE BY SUBROUTINE ( SUM2 ) 

(SUM3) 

(DELTA AND U3 BOTH DIFFERENT 
WHICH ASSUMES FROM ZERO ) 

( DELTA=0 BUT <J1+U2 AND U3 NOT 
EQUAL TO ZERO) 

(Ul, J2, AND DELTA ALL ZERO BUT 
U3 NOT ZERO) 

or 
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(SUM4 & SUM6) 

SUMMING WAS DONE BY SUBROUTINES (SUM4 & SUM6) 

(SUM5 & SUM6 ) 
(SUM5 & SUM6) 

(U3 = 0 BUT DELTA AND 01+02 NOT ZERO) 
ASSUMING (01,02, & 03 ALL ZERO BUT DELTA NOT 

ZERO) 

(03 AND DELTA ZERO BUT 01+02 NOT 
ZERO) 

(Ol, 02, 03 AND DELTA ALL ZERO) 



Following this is a statement of whether or not 
covergence test were met, and if not, a measure of the 



error introduced. 
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APPENDIX A 



PROGRAM LISTING 



A. IDENTIFICATION 



Title 

Co-op Id. 
Category- 
Programmer 
Date 



Electric Field Components by de Wette 
- Schacher Equations 
Zl-NPGS-LATSUM (F- 60 , 63 ) 

All Others 
David Dickmann 
February 1966 



B. PURPOSE 



This Fortram program computes electric fields per 
unit-dipoles for a lattice of similar dipoles. The 
program accepts lattices of any symmetry including 
triclinic and all higher symmetries, and computes field 
components in a cartesian reference frame, using the 
de Wette-Schacher equations. 

C. USAGE 



1. Calling sequence: This is a main program. See 

input formats. 

2. Input arguments: 

Al, A2, A3 are the lengths, in Angstroms, of 
the lattice unit cell basis vectors. 

XJ1, XJ2, XJ3 are fractions of the basis 
vectors Al, A2, A3 designating the location of 
the point at which the fields are to be evaluated. 

XADPH is the crystal angle between the A2 and 
A3 axes, in degrees and decimals. 

XBET is the crystal angle between the Al and 
A3 axes, in degrees and decimals. 

i 

XGAM is the crystal angle between the Al and 
A2 axes, in degrees and decimals. 
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3 . 
h * 

5 . 

6 . 



XDEL is the inclination— of the A3 axis from the 
perpendicular to the A1-A2 plane, in degrees and 
decimals. If unknown, enter as XDEL=90 and 
correct value is generated by program. 

AEX, AEY, AES are the absolute error limits 
for the values of the x, y, and z components of 
the electric fields. 

REX, REY, RES are the relative error limits 
for the corresponding electric field components. 

N is a summation cutoff value. If error limits 

p 

are not met after ^ 4N terms are summed, the 
program returns the field values at that time 
and a measure of the maximum error introduced 
by cutoff. 

Space required: 6446 cells (F-60), 4209 cells (E- 63 ) 

Temporary storage: None 

Alarms: None 

Error returns: 

a. Error in input parameters — Control returned 
to monitor with printout: 

PRO-CRAM STOPPED BY INPUT 

ERROR DESIGNATED BY K= _ _ 

Individual error codes are: 

K Error 



1 

2 

4 

10 

20 

40 

100 

200 



AU0. 

A2^ 0. 

A3^0. 

XJ1 <0. or XJ1 > 1. 

XJ2< 0. or X J2 > 1. 

XU3 <0. or XJ3 > 1. 
XALPHiO. or XALPH>90 
XBET^O. or XB ET > 90 . 
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1+00 XGAM^O. or XGAM^180. 

1000 XDEL<0. or XDEL>90. 

2000 N£0 

Multiple errors result in a K value which is 
the sum of the individual values above. 

b. N value too small to ensure good test for 

convergence of field sums. Control returned 
to monitor with self-explanatory printout. 

7. Error Stops: None except those associated with 

machine arithmetic errors. 

8. Input and output tape mountings: Normal 

9. Input and output formats 

a. Input is contained on four cards under the 
following format control: 

READ 1,A1,A2,A3,XvJ1,XvJ2,XJ3,XALFH, 

XB ET ,XGAM, XDEL ,N , A EX, REX, 

AEY , REY , AES ,RE2 

1 FORMAT (3F10.0/3F10.0AF10.0/I4,6E10.0) 

b. Output: Various 

(1) Error outputs as per 6 

(2) S.. - - the i-th component of the electric 

1J field due to a unit dipole oriented in 
the j-direction. Printed as an array 
under 1FIE19.10 with labels SUMXX, etc. 

(3) Hollerith prints indicating subroutines 
used and results of convergence tests. 

10. Selective stop and jump settings: Not applicable 

11. Timing: 

a. Compile time: 1+ min 20 sec 



21 



12. 



b. Run time: Average 1.4 sec per lattice, 

averaged over 336 combinations of lattices, 
error limits, and evaluation points. 

Accuracy: Relative error of S ; . ^min. .(REX,REY, 

RE2)*S.. J ’ J 

Absolute errdlr of Sj.-^minj :(AEX,AEY, 
AES) 

13. Equipment configuration: CDC 1604 with 

FORTRAN compiler or CO-OP- monitor system 
including F -63 

14. Library routines used: 

SQRTF 

SINF 

COSF 

ATANF 

ASINF 

ACOSF 

EXPF 

ABSF 

FIXFLOAT 

FLOAT FIX 

15* References: 

a. D. B. DICKMANN, A Computer Program 
for Evaluating Lattice Sums, Thesis, U.S. 
Naval Postgraduate School, Monterey, 1966 

b. F. W. de Wette and G. E. Schacher, Phys. 
Rev., 137, A78 (1965) 

D. METHOD OF ALGORITHM See references (a) and 
(b). 

E. FLOW CHART See reference (a). 
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APPENDIX B 



SUMMATION CONVERGENCE AND ERRORS 

The general lattice sums to be evaluated in the 
program are of the form 

S >J - Cl ^~u z hi > Ml )0! ('Un'Uz ) £ / er c * h t * 






e-^ h i > 

I -$3 (M, J 



■where C. is a constant, h. ,yU a ) is a function of 
and/A, which generally increases with increasing,/^ and 
and 0i is a phase factor of the form exp IrUTT i (CjyU, 
Because of the phase factors it is possible to have a 
distribution of significant terms over the plane such 
as shown in figure 8. 

To ensure inclusion of significant terms in such a 
distribution would require either very complex convergence 
tests or inclusion of many more terms then required. 

To avoid this, the sums were converted to semi- 
infinite sums over the first quadrant of th o plane. 
This led to the possible distributions shown in figure 9, 
where the degree of asymmetry about the axis 
depended on the ratio of the lengths of the a^ and a^ 
basis vectors and the angle between them, and the 
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m 





Figure 8: Possible distributions of significant terms of 

sum over pl ane « Shaded area contains 

all terms of absolute value larger than a given 
constant. 
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Figure 9 s Possible distributions of significant terms of 
converted sums 
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particular sum being evaluated, S-j. While these distributions 
still leave problems .in developing convergence tests, the 
problems are not as formidable as for the original equations. 

The solution used in this program proceeded in 
several steps. First the axial terms , 2, 

3, . . and^/Z^O, ^Z/^=l, 2,3- • • were summed together 

and the convergence test applied to the sum of each pair, 

> /J i //l=i> > until the tests were met for each 

sum Sjk, or N pairs were summed. 

The non-axial terms were then treated separately. 

First a skew factor, GIG = 2 * MAX (£?, costf+1, 

was generated. This number is such that the largest term 
of the sum (in absolute value) must occur at some point 
within the box Q&AJ/jMg, — J I <5 • This was tested against N 
to ensure N>GIG so that the sums would include the 
largest term. 

Summing was done by summing all terms /Jg = 1, 2, 

3 , . . . for ^1 = 1 and then for^/Z^=2, etc. The convergence 

test were threefold. First to ensure for given that 
assumed all values up to GIG and that this was done 
for all values of up to GIG. This ensured picking up 
the most significant term of this sum regardless of its 
magnitude. Second to ensure that for given /A-^=a. if sotne 
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term, ^2 — a > did not meet the convergence tests then 
summing would be done on the pencil /^=a+l , > 2, . . , 

a+1. Last, that summing would not exceed the terms 
lil= N, / 6/2 = i> /^i = j>^ 2 = N> regardless of other tests. 

The flow chart for this testing is shown in figure 10. 
Actual convergence zones for a cubic lattice and a triclinic 
lattice are shown in figures 11 and 12. 

After completing the program it was checked by 
running sums on 112 lattices (distortions of the basic cubic 
lattice) for which sums are known from geometric 
considerations. The error analysis for a relative error 
limit (REX, etc) of 10”\ 10”^, and 10”® and an absolute 
error limit (AEX, etc) equal to the relative limit is given 
in Table 2. 
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Figure 10a Flow chart of convergence tests common to 
SUM1-SUM6. Integer variable I refers to 
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Figure 10b Continuation of flow chart 




CONTINUE 

PROGRAM 
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Figure 11: Actual convergence zones plotted for unit 

cubic lattice using: library subroutine 
CONTOUR. 



il 
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Figure 12: Actual convergence zones plotted for a 

triclinic lattice using library subroutine 
CONTOUR . 
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Table 2: Error analysis of known sums the 

left-hand numbers are absolute errors, the right 
hand numbers are relative errors, the top numbers 
are maximum values, the bottom are average values, 
and the center integer is the number of sums 
evaluated. 



2a: 



Input relative error limit 
Input absolute error limit 
All errors X 10“° 




Lattice Subroutine Used 

Type SUM1 SUM2 SUM3 SUM4&6 SUM5&-6 



Triclinic 


.02 .1 
174 

.003 -005 


0 


0 


.1 .3 

156 

.008 .005 


0 


Monoclinic 


.03 .07 

12 

.004 .01 


.1 .2 

24 

.03 .01 


0 


.05 .09 

8 

.008 .01 


.04 .06 

16 

.01 .01 


Orthorhombic 


0 


1. .6 

18 

.08 .05 


0 


0 


.7 .4 

12 

.08 .05 


Tetragonal 


0 


.2 .1 

54 

.07 .06 


0 


0 


.09 .1 

36 

.003 .009 


Trigonal 


0 


0 


0 


0 


0 


Cubic 


0 


T5 7T 

18 

.07 .01 


0 


0 


.05 .08 
12 

.005 .007 
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2b: 



—6 

Input relative error limit = 10“ 
Input absolute error limit = 10“® 
All errors X 10~® 



Lattice Subroutine Used 

T yP e SUM1 SUM2 SUM3 SUM4&6 SUM5&6 



Triclinic 


.1 1. 

174 

.02 .03 


0 


0 


1 . 1. 
156 

.02 .01 


0 


Mono clinic 


.2 .5 

12 

.06 .1 


1. 1. 

24 

.1 .1 


0 


.3 1. 

8 

.05 . 2 


.2 .09 

16 

.03 .01 


Orthorhombic 


0 


4. 2. 

18 

.3 .2 


0 


0 


2. .9 

12 

.2 .1 


Tetragonal 


0 


2. 1. 

54 

.1 .1 


0 


0 


.8 1. 
36 

.04 .03 


Trigonal 


0 


0 


0 


0 


0 


Cubic 


0 


6 . .9 

18 

.4 .07 


0 


0 


.5 .9 

12 

.05 .1 
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2c: Input relative error limit = 10 

Input absolute error limit = 10 

All errors X lO - ^ 



-8 

-8 



Lattice 
Type 



SUM1 



Subroutine Used 
SUM2 



Triclinic 


.2 1. 

174 

.04 .03 


0 


0 


.9 1. 

156 

.04 .02 


0 


Mo no clinic 


.1 .4 

12 

.01 .04 


.9 1. 

24 

.05 .07 


0 


.1 .8 

8 

.02 .1 


.3 .08 
16 

.04 .007 


Orthorhombic 


0 


3. 1. 

18 

.2 .1 


0 


0 


2. 1. 

12 

.3 .1 


Tetragonal 


0 


3 . 2 . 

54 

.06 .05 


0 


0 


1. .9 

36 

.04 .03 


Trigonal 


0 


0 


0 


0 


0 


Cubic 


0 


7. 8. 

18 

.4 .5 


0 


0 


.4 -8 

12 

.1 .07 



67 



APPENDIX C 



INCOMPLETE GAMMA FUNCTION 

In evaluating the lattice sums when the evaluation 
point lines in the x-y plane an auxilliary convergence 
function is used to smooth and remove the simple pole of 
the function being summed. In the equations of Schacher 
and de Wette this introduces the following four terms into 
the sums: 

">». , Tth*) 

rc/i,nh x ) 
rm.nh*) 
r(% , Trh v 

where h=h.C^U, } Mz. ) > Mi ) sUz summing indices; and 

oo 

r(a,x) s/ t a -'e- r at the incomplete gamma function. 

In order to develop a usable program based on Schacher 
and de Wette's equation it was necessary to obtain a simple 
subroutine for evaluating these functions repeatedly as h 
changed with andylSg. Unfortunately incomplete gamma 
functions are not widely used and the only subroutines 
available from other sources were those designed to tabulate 
P (a, x) over a wide range of values of a and x. Since 
Pi a, x) is a transcedential whose "best" approximation is 
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an assymptotic series, tabulating programs are both 
cumbersome and slow. 

Three possibilities were considered for this program. 
First, adapt one of the existing programs for limited use. 

This was tried but even after abstracting the desired 
parts from the best of these programs we were left with 
a double-precision program with an evaluation time, per h 
value, of about 0.5 sec. Since we contemplated hundreds 
of thousands of evaluations per lattice we ultimately 
discarded this approach. 

The second possibility considered was simply storing 
tabulated values over the range h=0-100 and interpolating 
as necessary. However to obtain program accuracies of 
±10"9 f or Sjj would have required storing 4000 or more 
values. This, added to the large program storage already 
generated would have overflowed the internal storage of 
all but the largest computers available. Tape or disc 
storage would again increase run time too much so this 
approach was also discarded. 

The final possibility was adopted, that of formulating 
new approximations to the incomplete gamma function. 
Ultimately a set of power series approximations was developed, 
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7 

using the Tau Method of Lanczos, to span the range 
77 h* =3-00. Another, long-known series, well convergent 
for small values of x, was used in the range 7f/)* =0-3 



The expansions are 








r('/2>xh 




3 < x 



where the C^i s were generated using the Tau Method. The 
other functions, , etc were obtained from rC/z^x) 

by using the recursion relation r (0 + / t X) -ClFCCI)#) +X a 0~* * 
The Tau Method is well adopted to exploring several 
possibilities of approximating series on the computer. In 
this instance we examined a total of 240 possibilities, 
combining the four a values, -1, J, £ with 3 possible 
ranges of approximation, j-<X , !<X } 3<X and in each case 
examining series of 2 , 3> 4 • ••, 21 terms. In general, 
accuracy increased with increasing a value, increasing end 
point b, b< x, and increasing number of terms in the series. 
We used a=j£ vice j| or £ because it matched the a value of 

7 

Lanczos, C., Applied Analysis , 464-485, (Prentice 
Hall, Englewood Cliffs, N.J., 1961) • 
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the expansion for 0<x^3» -A. 17 term series was used 



since it exhausted machine accuracy of the computer used 
*“11 

(±10 ), and any higher series would add no accuracy. 

Having settled on an approximation it was checked 
against a double precision routine known to be accurate to 
±10”'*''*'. The graph of the error of the approximation is 
shown in figure 13 . 



71 




Figure 13! Error curve for Incomplete Gamma Function 
approximation. Abscissa is x, ordinate is 
r* {'/ZiXl-r'C/x,*) where/** is approximation and 
P is actual value of Incomplete Gamma Function. 
Shaded area represents an oscillation of error 
of such a high fi'equency that it could not be 
faithfully represented on this scale. 
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APPENDIX D 



APPLICATION TO BENZENE 

During the preparation of this program R.M. Hexter 
of the Mellon Institute suggested an application to a long- 
standing problem in molecular spectroscopy. This work has 
been completed and is detailed below as an example of applying 
the program to an actual problem. 

In calculating the electronic spectrum of benzene it 
is necessary to compute the transition moments of the free 
molecules. These moments depend on the dipole-dipole sums 
of the crystal. Various workers computing these sums by 

various methods have obtained varying values for the sums. 

8 

Craig and Walsh have set out the crystalline structure 
assumed and tabulated the values obtained by two summation 
procedures. This program was used to check these values 
with results tabulated below. 

The benzene crystal model consists of an orthorhombic 
unit cell with four molecular orientations at the corners 

and faces. If 'we superimpose cartesian reference frames 

£ 

on each molecule we can identify our molecules as follows: 

®D. Craig and tJ. Walsh, <J . Chem Soc 1958, 1613 
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Molecular Unit Cell 


Crystal 


Direction Cosines 


Type Location 


Axis 


X 


y 


z 


I 0,0,0 


a 


.6365 


-.2992 


.7109 




b 


.2000 


.9541 


.2229 




c 


-.7450 


0 


.6671 


11 i-°’i 


a 


-.6365 


.2992 


-.7109 




b 


-.2000 


-.9541 


-.2229 




c 


-.7450 


0 


.6671 


111 °>i>i 


a 


-.6365 


.2992 


-.7109 




b 


.2000 


.9541 


-.2229 




c 


• 7450 


0 


.6671 


IV i>i>° 


a 


.6365 


-.2992 


.7109 




b 


-.2000 


-.9541 


-.2229 




c 


.7450 


0 


-.6671 


Where (a,b,c) forms the cartesian frame 


of the 


unit cell. 



The unit cell parameters are (at -3° C) 

a x = 7.46O A or 

ag = 9 • 666 A. 

a 3 = 7.034 A 

The problem then is to compute the interaction energy 
of a unit dipole in the i- direction in molecule I due to a 
lattice of unit dipoles in the j-direction in molecules of 
type V , where i, j = x, y, z; V = I, II, III, IV. 

In comparing figures obtained by this program with 
those tabulated by Craig and Walsh two difficulties arise. 
First dipole-dipole energies from this program work out in 
recipricol angstroms per coulomb squared while Craig and 
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Walsh tabulate in recipricol centimeters (wave number) per 
electronic charge squared. This requires multiplying the 
results of this program by 1.161 X 10^ to be comparable. 

Second Craig and Walsh assume an- infinite* sphSH'Slfl 
crystal while this program assumes an infinite slab-shaped 
crystal. Thus depolarization factors must be applied. The 
two systems are: 

Sperical Slab-shaped 



r- 4TT 

3 Vo/ 




r— 4-Tf 

£yy 3Vol 




- sirr 

lZ 2Z 3VcT 


r- 47T' 

£-22 Vo/ 



Thus, to be comparable we must add - (l+7r/vol) = 
0.0165170 to zT 2 !;' ! terms generated by this program and 
subtract - (l+77'/vol) = 0.0082585 to xx and xy terms 
generated by this program. (Where vol = volume of unit 
cell in A.) 

GENERAL PROCEDURE 
With these considerations in mind the computation 
becomes straight-forward if tedious. First we will compute 
electric field components at molecule I along the crystal 
axes a, b, c due to an assumed set of unit dipoles along 
the a, b, and c axes at molecular positions I, II, III, and 
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IV. This is done by calculating the field at the appropriate 
position of molecule I (given by XJI, XJ2, X.J3) within 
each of the I, II, III, and IV sublattice unit cells. Then 
we correct the depolarization factors and apply the conversion 
factor. Then, using the direction cosines of the type I 
molecule for each molecular axis x, y, and z (with respect 
to the a, b, c axes) we can get the x, y, and z components 
of field at I due to a, b, and c dipoles at all molecular 
positions. Then by multiplying by the direction cosines for 
each molecule respectively we get energies of the x, y, and 
z dipoles of molecule I due to x, y, and z dipoles at each 
molecular position. These we compare with similarly 
tabulated results of Craig and Walsh. 

DETAILS OF COMPUTATION 

To get the a, b, and c field components at molecule 
I we enter the program with: 

A1 = 7.460 XALPH = 90. 

A2 = 9.666 XBET = 90. 

A 3 = 7.034 XGAM = 90. 

XDEL = 0. 

AEX=REX=AEY=REY=AE2=RE2=10 -7 

t 

N=20 

and for various molecules: 
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X J1 
XJ2 
XJ3 



I II III IV 

o ±04 

o o i 4 

o 4 4 o 



Since the crystal is orthorhombic we get S.j=0, i^j. 



The three nonzero outputs SUMXX, SUMYY, and SUM22 
here correspond to the a, b, and c field components. The 
computed Sjj values, depolarization correction, conversion 
factor, and converted values are listed below: 



Fields at 


i-f ield 








I due to 


component 








molecules 


due to j- 


Computed 




Corrected 


of type 


dipoles 


Values 


Depolarization 


Values 


I 


aa 


.0104233 


-.0082585 


.0021648 




bb 


.0011720 


-.0082585 


.0070865 




cc 


.0115953 


.0165170 


.0049217 


II 


aa 


.0240220 


-.0082585 


.0157635 




bb 


-.0187985 


-.0082585 


-.0270570 




cc 


-.0052235 


.0165170 


.0112935 


III 


aa 


-.0055470 


-.0082585 


-10138055 




bb 


.0231906 


-.0082585 


.0149321 




cc 


-.0176436 


.0165170 


-.0011266 


IV 


aa 


.0080050 


-.0082585 


-.0002535 




bb 


.0209184 


-.0082585 


.0126599 




cc 


-.0289234 


.0165170 


-.0124064 
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Corrected 


Conversion 


Converted 


Molecule 


Axes 


Value 


Factor 


Values 


I 


aa 


.002161+8 


1. 161X10 5 


251.333 




bb 


-.0070865 






-822.71+2 




cc 


.001+9217 






571.1+09 


II 


aa 


.0157635 






1830.11+2 




bb 


-.0270570 






-3H+I.3I7 




cc 


.0112935 






1311.175 


III 


aa 


-.0138055 






-1602.818 




bb 


. 011+9321 






1733. 616 




cc 


-.0011266 






-130.798 


IV 


aa 


-.0002535 






-29.1+31 




bb 


.0126599 


1 


1 _ 


11+69.811+ 




cc 


-.0121+061+ 


1. 161X10^ 


-11+1+0.383 



Now we need the x, y, and z components of these 



fields at molecule I; I-I (x, m), I-I (y, m) , I-l (z, m), 
etc (m=a,b,c). Here we need the molecule I direction 
cosines for x, y, or z with respect to a, b, and c* Two 
examples suffice: 



x - field components 



Fields at I 



due to dipoles 
at 


Unit Call 
Component 


Value 


Dir. Cos 


Molecular 

Component 


I 


aa 


251.333 


.6365 


159.9734 




bb 


-822.742 


.2000 


-164-5484 




cc 


571.1+09 


-.71+50 


-425-6997 



II ' 

t t 

i i 
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y - field components 



Fields at I 

due to dipoles Unit Call Molecular 

at Component Value Dir Cos Component 

i 

t 

i 



aa 


1830 . 142 


-.2992 


-547. 


5784 


bb 


-311+1.317 


• 9541 


-2997 


.1305 


cc 


1311.175 


0 


0 





III 

i 

i 

From these we can make up the energies by multiplying 
the molecular components by the direction cosines of the 
source dipole at the molecular position of interest. 

Since a we have tabulated the molecular components 
I-I (x, m) and I-II (y, m), if we want the energy of an 
x-dipole in molecule I due to a lattice of y-dipoles at 
position I, I-I (x,y), we need only multiply the molecular 
components I-I (x, m) by the y direction cosines of type 
I and sum. We show the computation of I-I (x,y) and 
I-II (y,z) below: 
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Field Direction Energy 

Component Axis Value Cosine Contribution 

Iy dipoles 

I-I (x, ) aa 159.9734 -.2992 -47.8640 

bb -164.5484 .9541 -156.9956 

cc -425.6997 0 0 

I-I (x,y) energy -204*8596 

II z dipoles 

I-II (y, ) aa -547-5784 -.7109 389.2734 

bb-2997.1305 -.2229 668.0603 

cc 0 .6671 0 

I-II (y,z) energy 1057.3337 

These then are the numbers to be compared with 

those tabulated by Craig and Walsh. Their table contains 

q 

numbers obtained both from the Ewald-Kornf eld 7 method 
applied to an infinite sphere and from a direct molecule- 
molecule summation over a sphere of radius 20 A. The 
comparison is shown below. Results (a) are from this 
program, (b) from the Ewald-Kornf eld method, and (c) 
from direct molecule-molecule summation. 



9 

Kittel, Introduction to Solid State Physics, 2nd Ed, 
347, (Cohn Wiley and Sons, Inc., New York, 1954) 
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Dipole 

Axes Results I, I 



X, X 


(a) 


-386.06 




(b) 


-386 




(c) 


-415 


y» y 


(a) 


726.45 




(b) 


727 




(c) 


790 


z , z 


(a) 


-340.47 




(b) 


-341 




(c) 


-375 


X, y 


(a) 


204-86 




(b) 


205 




(c) 


233 


X, z 


(a) 


206.94 




(b) 


207 




(c) 


181 


y> z 


(a) 


228.43 




(b) 


228 




(c) 


260 



Molecule 
I, II 


Sites 
I, III 


I, IV 


-111.94 

-112 

-92 


-791.29 

-791 

-769 


-728.72 

-729 

-771 


-2695-73 

-2679 

-2629 


-1721. 61 

-1723 

-1713 


1340-62 

1341 

1405 


185-49 

185 

220 


-954-10 

-955 

-950 


-553.20 

-553 

-587 


-947-96 

-947 

-934 


-25-57 

-25 

-8 


274-87 

275 

289 


1339-72 

1340 

1306 


-737-54 

-738 

-818 


794-70 

795 

833 


-1057-33 

-1058 

-1041 


-27-76 

-27 

-7 


306.32 

307 

323 
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